/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
#include "EiccFdt.h"

#include "FairGeoInterface.h"
#include "FairGeoLoader.h"
#include "FairGeoNode.h"
#include "FairGeoVolume.h"
#include "FairRootManager.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairVolume.h"
#include "EiccDetectorList.h"
#include "EiccStack.h"
#include "EiccFdtGeo.h"
#include "EiccFdtGeoPar.h"
#include "EiccFdtPoint.h"

#include <TClonesArray.h>
#include <TGeoBBox.h>
#include <TGeoCompositeShape.h>
#include <TGeoManager.h>
#include <TGeoMaterial.h>
#include <TGeoMedium.h>
#include <TGeoTube.h>
#include <TGDMLParse.h>
#include <TVirtualMC.h>
#include <iostream>
using std::cout;
using std::endl;

EiccFdt::EiccFdt()
    : FairDetector("EiccFdt", kTRUE, kEiccFdt)
    , fTrackID(-1)
    , fEventID(-1)
    , fVolumeID(-1)
    , fPos()
    , fMom()
    , fTime(-1.)
    , fLength(-1.)
    , fELoss(-1)
    , fEiccFdtPointCollection(new TClonesArray("EiccFdtPoint"))
{}

EiccFdt::EiccFdt(const char* name, Bool_t active)
    : FairDetector(name, active, kEiccFdt)
    , fTrackID(-1)
    , fEventID(-1)
    , fVolumeID(-1)
    , fPos()
    , fMom()
    , fTime(-1.)
    , fLength(-1.)
    , fELoss(-1)
    , fEiccFdtPointCollection(new TClonesArray("EiccFdtPoint"))
{}

EiccFdt::EiccFdt(const EiccFdt& right)
    : FairDetector(right)
    , fTrackID(-1)
    , fEventID(-1)
    , fVolumeID(-1)
    , fPos()
    , fMom()
    , fTime(-1.)
    , fLength(-1.)
    , fELoss(-1)
    , fEiccFdtPointCollection(new TClonesArray("EiccFdtPoint"))
{}

EiccFdt::~EiccFdt()
{
    if (fEiccFdtPointCollection) {
        fEiccFdtPointCollection->Delete();
        delete fEiccFdtPointCollection;
    }
}

void EiccFdt::Initialize()
{
    /**
     * WORKAROUND needed for Geant4 in MT mode
     * Call AddSensitiveVolume for sensitive volumes in order to fill
     * thread-local FairModule::svList.
     */
    DefineSensitiveVolumes();

    FairDetector::Initialize();
    FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
    EiccFdtGeoPar* par = (EiccFdtGeoPar*)(rtdb->getContainer("EiccFdtGeoPar"));
}

Bool_t EiccFdt::ProcessHits(FairVolume* vol)
{
    /** This method is called from the MC stepping */

    // Set parameters at entrance of volume. Reset ELoss.
    if (TVirtualMC::GetMC()->IsTrackEntering()) {
        fELoss = 0.;
        fTime = TVirtualMC::GetMC()->TrackTime() * 1.0e09;
        fLength = TVirtualMC::GetMC()->TrackLength();
        TVirtualMC::GetMC()->TrackPosition(fPos);
        TVirtualMC::GetMC()->TrackMomentum(fMom);
    }

    // Sum energy loss for all steps in the active volume
    fELoss += TVirtualMC::GetMC()->Edep();

    // Create EiccFdtPoint at exit of active volume
    if ( (TVirtualMC::GetMC()->IsTrackExiting() || TVirtualMC::GetMC()->IsTrackStop()
        || TVirtualMC::GetMC()->IsTrackDisappeared())&& TVirtualMC::GetMC()->TrackCharge()!=0 ) {
    //if ( (TVirtualMC::GetMC()->IsTrackExiting() || TVirtualMC::GetMC()->IsTrackStop()
    //    || TVirtualMC::GetMC()->IsTrackDisappeared()) ) {
        fTrackID = TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber();
        fVolumeID = vol->getMCid();
	int copyNo= -1;
	TVirtualMC::GetMC()->CurrentVolID(copyNo);
	//std::cout<<"=====================================name:  "<<TVirtualMC::GetMC()->CurrentVolName()<<"    "<<std::endl;
	int PID  = TVirtualMC::GetMC()->TrackPid();
	fTime = PID;

	//std::cout<<"XYZ are "<<fPos.X()<<", "<<fPos.Y()<<", "<<fPos.Z()<<",  R is "<<sqrt(fPos.X()*fPos.X()+fPos.Y()*fPos.Y())<<", ID is "<<fVolumeID<<std::endl;
	fEventID = gMC->CurrentEvent(); //yutie 0407
        if (fELoss == 0.) {
            return kFALSE;
        }

	//std::cout<<"XYZ are "<<fPos.X()<<", "<<fPos.Y()<<", "<<fPos.Z()<<std::endl;
        AddHit(fTrackID,
               copyNo,
               TVector3(fPos.X(), fPos.Y(), fPos.Z()),
               TVector3(fMom.Px(), fMom.Py(), fMom.Pz()),
               fTime,
               fLength,
               fELoss,
	       fEventID); 

        // Increment number of EiccFdt det points in TParticle
        EiccStack* stack = (EiccStack*)TVirtualMC::GetMC()->GetStack();
        stack->AddPoint(kEiccFdt);
    }

    return kTRUE;
}

void EiccFdt::EndOfEvent()
{

    LOG(info) << "EiccFdt: " << fEiccFdtPointCollection->GetEntriesFast() << " points registered in this event";

    fEiccFdtPointCollection->Clear();
}

void EiccFdt::Register()
{

    /** This will create a branch in the output tree called
      EiccFdtPoint, setting the last parameter to kFALSE means:
      this collection will not be written to the file, it will exist
      only during the simulation.
  */

    if (!gMC->IsMT()) {
        FairRootManager::Instance()->Register("EiccFdtPoint", "EiccFdt", fEiccFdtPointCollection, kTRUE);
    } else {
        FairRootManager::Instance()->RegisterAny("EiccFdtPoint", fEiccFdtPointCollection, kTRUE);
    }
}

TClonesArray* EiccFdt::GetCollection(Int_t iColl) const
{
    if (iColl == 0) {
        return fEiccFdtPointCollection;
    } else {
        return NULL;
    }
}

void EiccFdt::Reset() { fEiccFdtPointCollection->Clear(); }

void EiccFdt::ConstructGeometry()
{

    TString fileName=GetGeometryFileName();
    if(fileName.EndsWith(".root")) {
	    ConstructRootGeometry();
    /*} else if(fileName.EndsWith(".gdml")) {
	    ConstructGDMLGeometry();*/
    } else {
	    ConstructAsciiGeometry();
    }

}

/*
void EiccFdt::ConstructGDMLGeometry(){

    //Open GDML file in EiccRoot/geometry   "genfitGeom_AllSi_v3.gdml"
    TString fileName=GetGeometryFileName().Data();
    TGeoVolume* top = gGeoManager->GetTopVolume();
    TGDMLParse parser;
    TGeoVolume *dsk_with_pipe = parser.GDMLReadFile(fileName);
    if(dsk_with_pipe){
	    top->AddNode(dsk_with_pipe, 0);
    }
}
*/

void EiccFdt::ConstructAsciiGeometry(){

    TGeoVolume* top = gGeoManager->GetTopVolume();
    TGeoMedium* Si = gGeoManager->GetMedium("silicon");
    TGeoMedium* Carbon = gGeoManager->GetMedium("carbon");
    TGeoMedium* Aluminium = gGeoManager->GetMedium("Aluminium"); //26.98 13. 2.7
    TGeoMedium* Air = gGeoManager->GetMedium("air");
    TGeoMedium* PWO = gGeoManager->GetMedium("PWO");

    if (Si == 0) {
	    std::cout<<"no si???"<<endl;
        TGeoMaterial* matSi = new TGeoMaterial("Si", 28.0855, 14, 2.33);
        Si = new TGeoMedium("Si", 2, matSi);
    }
    if (Carbon == 0) {
	    cout<<"no carbon???"<<endl;
        TGeoMaterial* matCarbon = new TGeoMaterial("C", 12.011, 6.0, 2.265);
        Carbon = new TGeoMedium("C", 3, matCarbon);
    }
    if (Aluminium == 0) {
	    cout<<"no aluminum???"<<endl;
        TGeoMaterial* matAlum = new TGeoMaterial("Aluminium", 26.98, 13., 2.7);
        Aluminium = new TGeoMedium("Aluminium", 4, matAlum);
    }



    double thickness = 0.005;//0.005cm  =  0.00005m = 50 um
    double thickness_alum = 0.02; //set alum thickness to be 4 time of the pixel chip.

    /////
    //
    /*
    for(int i = 0; i < 10; i++){
	TGeoVolume* temp_dipole2 = gGeoManager->MakeBox(Form("temp_%d", i+1), Si, 10, 10, thickness/2.);
	TGeoRotation r_temp;
        r_temp.SetAngles(0,0,0);
        r_temp.RotateY(1.0*(-1*0.065)*57.2958);
        TGeoTranslation t_temp(-53+i*15*sin(-0.065), 0, 850.+i*15);
        TGeoCombiTrans c_temp(t_temp, r_temp);
        TGeoHMatrix* h_temp = new TGeoHMatrix(c_temp);
        top->AddNode(temp_dipole2, i+1, h_temp);
        AddSensitiveVolume(temp_dipole2);
    }
*/

    double RomanPotHalfX = 18; //18;
    double RomanPotHalfY = 6.; //6;
    double RomanPotHoleHalfX = 3.8;//6.0;
    double RomanPotHoleHalfY = 1.7;//3.0;
    double RomanPotAngle = -0.050-0.015-0.070;

    const int NRomanPot = 3;
    double RomanPotZ[NRomanPot] = {1000. + 15.*cos(RomanPotAngle), 1000. + 25.*cos(RomanPotAngle), 1000. + 35.*cos(RomanPotAngle)};
    double RomanPotX[NRomanPot] = {-63.2345 + 15.*sin(RomanPotAngle), -63.2345 + 25.*sin(RomanPotAngle), -63.2345 + 35.*sin(RomanPotAngle)};

    std::cout<<"fdt disk info:"<<std::endl;

    TGeoTranslation *RomanpotM1 = new TGeoTranslation(0, 0, 0);
    RomanpotM1->SetName("m1_RomanPot");
    RomanpotM1->RegisterYourself();

    for (int i=0; i<NRomanPot; i++){
	std::cout<<RomanPotZ[i]<<" "<<RomanPotX[i]<<std::endl;
        TGeoVolume* RomanPotVol1 = gGeoManager->MakeBox(Form("Roman_pot_Vol1_%d", i+1), Si, RomanPotHalfX, RomanPotHalfY, thickness/2.);
        TGeoVolume* RomanPotVol2 = gGeoManager->MakeBox(Form("Roman_pot_Vol2_%d", i+1), Si, RomanPotHoleHalfX, RomanPotHoleHalfY, thickness/2.*1.0001);

        TGeoCompositeShape *shape_RomanPotVol = new TGeoCompositeShape(Form("shape_RomanPotVol%d", i+1), Form("Roman_pot_Vol1_%d - Roman_pot_Vol2_%d:m1_RomanPot", i+1, i+1));
        TGeoVolume* RomanPotVol = new TGeoVolume(Form("Det_RomanPot_%d", i+1), shape_RomanPotVol, Si);
        TGeoRotation r_RomanPotTube;
        r_RomanPotTube.SetAngles(0,0,0);
        r_RomanPotTube.RotateY(1.0*RomanPotAngle*57.2958);
        TGeoTranslation t_RomanPotTube(RomanPotX[i], 0, RomanPotZ[i]);
        TGeoCombiTrans c_RomanPotTube(t_RomanPotTube, r_RomanPotTube);
        TGeoHMatrix* h_RomanPotTube = new TGeoHMatrix(c_RomanPotTube);
        top->AddNode(RomanPotVol, i+1, h_RomanPotTube);
        AddSensitiveVolume(RomanPotVol);

        TGeoVolume* RomanPotVolAx1 = gGeoManager->MakeBox(Form("Appendix1_roman_pot_%d", i+1), Aluminium, RomanPotHalfX, RomanPotHalfY, thickness_alum/2.*0.999);
	TGeoVolume* RomanPotVolAx2 = gGeoManager->MakeBox(Form("Appendix2_roman_pot_%d", i+1), Aluminium, RomanPotHoleHalfX, RomanPotHoleHalfY, thickness_alum/2.);
        TGeoCompositeShape *shape_RomanPotAxVol = new TGeoCompositeShape(Form("shape_RomanPotAxVol%d", i+1), Form("Appendix1_roman_pot_%d - Appendix2_roman_pot_%d:m1_RomanPot", i+1, i+1));
        TGeoVolume* RomanPotVolAx = new TGeoVolume(Form("Appendix1_roman_pot_%d", i+1), shape_RomanPotAxVol, Aluminium);

        TGeoTranslation t_RomanPotTubeAx(RomanPotX[i] + (thickness/2. + thickness_alum/2.)*sin(RomanPotAngle), 0, RomanPotZ[i] + (thickness/2. + thickness_alum/2.)*cos(RomanPotAngle));
        TGeoCombiTrans c_RomanPotTubeAx(t_RomanPotTubeAx, r_RomanPotTube);
        TGeoHMatrix* h_RomanPotTubeAx = new TGeoHMatrix(c_RomanPotTubeAx);
        top->AddNode(RomanPotVolAx, i+1, h_RomanPotTubeAx);
    }


    return;

    //calc det. info for the fdt disks
    double z_fdt = 1608; //z position = 12.0m
    double l_fdt = 80.;//200; //length = 1.5 m
    double angle_cover_fdt = 0.02;//0.050*1.55;

    double B_field = 3.4;//2.69; //Tesla
    double R_bending = 20*10./3./B_field * 100; //=26.455meter

    double angle_cover_dist = 0.019; //edt from 16 mrad. we set fdt maximum to 19 mrad.

    double rmax_ele_beam = 5.8; //5.5 + beampipe 1cm ?
    double rmax_ion_beam = 6; //5cm?
    double rmin_ele_beam = 5.5;
    double rmin_ion_beam = 5.5;
    double fdt_rmin1 = 0, fdt_rmax1 = (z_fdt-0.5*l_fdt) * atan(angle_cover_fdt);
    double fdt_rmin2 = 0, fdt_rmax2 = (z_fdt+0.5*l_fdt) * atan(angle_cover_fdt);

    //hard code of location of the dipole center: x1789.28 z=1808.93  R=1960.78
    //double x_center_dipole = 1789.28, z_center_dipole = 1808.93, R_dipole = R_bending;
    double x_center_dipole = 1801.65, z_center_dipole = 1861.97, R_dipole = R_bending;
    int N_disk = 4;
    double phi_fdt_coverage = 0.04; //phi coverage of FDT = 40 mrad
    double phi_fdt_edge = 0.005;//5mrad not used

    TGeoVolume* fdt_container = gGeoManager->MakeCone("fdt_endcap",Air, l_fdt/2., fdt_rmin1, fdt_rmax1, fdt_rmin2, fdt_rmax2);

    TGeoVolume *dipole_cole = gGeoManager->MakeCone("dipole_cole",Aluminium, l_fdt/2., fdt_rmax1*0.85, fdt_rmax1*0.9, fdt_rmax2*0.85, fdt_rmax2*0.9);
    TGeoRotation r_dipole;
    r_dipole.SetAngles(0, 0, 0);
    TGeoTranslation t_dipole(0., 0, 0);
    TGeoCombiTrans c_dipole(t_dipole, r_dipole);
    TGeoHMatrix* h_dipole = new TGeoHMatrix(c_dipole);
    //fdt_container->AddNode(dipole_cole, 1, h_dipole);
    dipole_cole->SetLineColor(kGreen+1);

    int det_color = kOrange+1;
    int app_color = kSpring-6;

    double expand_for_plotting = 0.9999;

    TGeoVolume* det_fdt_dsk[4]; //N_disk = 4
    TGeoVolume* Appendix_fdt_dsk[4];//N_disk = 4


    double phi_fdt_disk_i, x_disk_i, z_disk_i, rmin_disk_i, rmax_disk_i;
    for(int i = 0; i < N_disk; i++){
            phi_fdt_disk_i = -0.130 + phi_fdt_edge/2. + i*(phi_fdt_coverage-phi_fdt_edge)/(N_disk-1);
            x_disk_i = x_center_dipole - R_dipole*cos(phi_fdt_disk_i);
            z_disk_i = z_center_dipole + R_dipole*sin(phi_fdt_disk_i);
	    std::cout<<x_disk_i<<" "<<z_disk_i<<std::endl;
	    rmin_disk_i = rmax_ion_beam;
	    rmax_disk_i = z_disk_i*atan(angle_cover_dist);

	    det_fdt_dsk[i] = gGeoManager->MakeTubs(Form("Det_fdt_dsk_p%d", i+1), Si, rmin_disk_i, rmax_disk_i, thickness/2, 0, 360);

	    TGeoRotation r_fdt_i;
	    r_fdt_i.SetAngles(0, 0, 0);
	    TGeoTranslation t_fdt_i(x_disk_i, 0, z_disk_i);
	    TGeoCombiTrans c_fdt_i(t_fdt_i, r_fdt_i);
	    TGeoHMatrix* h_fdt_i = new TGeoHMatrix(c_fdt_i);
	    //fdt_container->AddNode(det_fdt_dsk[i], 1, h_fdt_i);
	    //top->AddNode(det_fdt_dsk[i], 1+i, h_fdt_i);
	    det_fdt_dsk[i]->SetLineColor(det_color);
	    //AddSensitiveVolume(det_fdt_dsk[i]);

	    Appendix_fdt_dsk[i] = gGeoManager->MakeTubs(Form("Appendix_fdt_dsk_p%d", i+1), Aluminium, rmin_disk_i, rmax_disk_i, thickness_alum/2*0.999, 0, 360);

	    TGeoTranslation t_fdt_appendix(x_disk_i, 0, z_disk_i+thickness_alum/2);
	    TGeoCombiTrans c_fdt_appendix(t_fdt_appendix, r_fdt_i);
	    TGeoHMatrix* h_fdt_appendix = new TGeoHMatrix(c_fdt_appendix);
	    //fdt_container->AddNode(Appendix_fdt_dsk[i], 1, h_fdt_appendix);
	    //top->AddNode(Appendix_fdt_dsk[i], 1+i, h_fdt_appendix);
	    //Appendix_fdt_dsk[i]->SetLineColor(app_color);

    }
/*
    TGeoRotation r_fdt;
    r_fdt.SetAngles(0,0,0);
    TGeoTranslation t_fdt(0, 0, 0);
    TGeoCombiTrans c_fdt(t_fdt, r_fdt);
    TGeoHMatrix* h_fdt = new TGeoHMatrix(c_fdt);
    top->AddNode(fdt_container, 1, h_fdt);
*/
}


EiccFdtPoint* EiccFdt::AddHit(Int_t trackID,
                                      Int_t detID,
                                      TVector3 pos,
                                      TVector3 mom,
                                      Double_t time,
                                      Double_t length,
                                      Double_t eLoss,
				      Int_t eventID)
{
    TClonesArray& clref = *fEiccFdtPointCollection;
    Int_t size = clref.GetEntriesFast();
    return new (clref[size]) EiccFdtPoint(trackID, detID, pos, mom, time, length, eLoss, eventID);
}

FairModule* EiccFdt::CloneModule() const { return new EiccFdt(*this); }

void EiccFdt::DefineSensitiveVolumes()
{
    TObjArray* volumes = gGeoManager->GetListOfVolumes();
    TIter next(volumes);
    TGeoVolume* volume;
    while ((volume = static_cast<TGeoVolume*>(next()))) {
        if (IsSensitive(volume->GetName())) {
            LOG(debug2) << "Sensitive Volume " << volume->GetName();
	    std::cout<<"Sensitive Volume " << volume->GetName()<<std::endl;
            AddSensitiveVolume(volume);
        }
    }
}

Bool_t EiccFdt::IsSensitive(const std::string& name)
{
    if (name.find("Det_fdt") != std::string::npos) {
        return kTRUE;
    }
    return kFALSE;
}

ClassImp(EiccFdt);
